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CHRISTINE:  A  MULTIFREQUENCY  PARAMETRIC  SIMULATION 
CODE  FOR  TRAVELING  WAVE  TUBE  AMPLIFIERS 


1.  INTRODUCTION 

In  many  applications  of  traveling  wave  tube  amplifiers  (TWTs)  in  the  defense  arena,  including 
electronic  warfare  (EW),  radar,  and  communications,  TWTs  are  required  to  amplify  multiple 
frequencies  simultaneously.  In  these  applications,  nonlinearities  in  the  amplification  process  cause 
intermodulation  and  amplitude  and  phase  cross-modulation.  When  the  number  of  carriers  is  large, 
this  phenomenon  becomes  even  more  complex.  Existing  tubes  avoid  these  parasitic  effects  by 
combining  such  techniques  as  feed-forward  compensation  or  operation  in  the  linear  regime  with 
depressed  energy  collection  of  the  electron  beam.  Associated  with  these  techniques  is  the  cost  of 
increasing  the  complexity  of  the  tube  or  the  weight  of  the  power  supply.  Some  EW  applications 
require  intermodulation  products  at  levels  no  greater  than  -30  to  -40  dBc  and  phase  stability  less  than 
3°/dB.  Communication  applications  can  require  intermodulation  levels  as  low  as  -70  dBc.  Moreover, 
the  requirement  to  further  increase  the  efficiency,  frequency,  and  power  level  of  TWTs  makes 
attaining  these  levels  of  intermodulation  and  cross-modulation  even  more  difficult.  Thus,  interest  has 
been  refocused  on  designing  TWTs  that  minimize  intermodulation  and  cross-modulation  distortion. 

This  report  describes  the  workings  of  a  parametric  computer  code,  CHRISTINE,  that  models  the 
operation  of  TWTs  in  the  presence  of  multiple  frequency  signals.  The  report  provides  examples  of 
how  CHRISTINE  can  be  used  and  serves  as  a  user’s  guide  to  the  model.  The  basic  model  is  a 
variation  of  that  which  has  been  in  use  in  the  field  for  some  time  [1-4].  We  have  attempted  to  keep 
the  notation  as  similar  as  possible  to  that  of  the  basic  equations  of  classical  physics.  The  parametric 
description  of  the  complex  interaction  between  the  particles  and  fields  in  a  TWT  is  not  all  inclusive  as 
is,  say,  that  of  a  particle  in  cell  code.  However,  CHRISTINE  can  quickly  calculate  the  output  spectrum 
of  a  device  with  many  frequencies  present.  As  such,  it  is  useful  in  designing  devices  where  many 
simulations  are  required  to  cover  a  large  parameter  space. 

Section  2  of  this  report  describes  the  physical  model,  which  forms  the  basis  for  the  code.  The 
discussion  is  general  enough  to  provide  a  user  of  the  code  a  good  understanding  of  which  processes 
are  being  simulated.  Section  3  describes  some  of  the  numerical  methods  used  in  solving  the  basic 
equations  and  discusses  the  organization  of  the  code  and  the  functions  of  the  various  subroutines. 
Section  4  provides  a  set  of  examples  that  illustrate  the  workings  of  the  code,  the  input  parameter 
namelists  and  files,  and  the  output  files. 

2.  PHYSICAL  MODEL 

2.1  The  Electromagnetic  Fields  of  the  Structure  and  Their  Frequencies 

We  begin  the  description  of  the  physical  model  employed  in  CHRISTINE  with  the  representation 
of  the  assumed  fields  of  the  cold  stmcture.  Ultimately,  in  the  type  of  model  presented  here,  only  two 
frequency  dependent  quantities — the  phase  velocity  and  the  coupling  impedance — are  needed  to 
specify  the  parameters  of  the  cold  structure  fields.  We  begin,  however,  at  a  more  basic  level  to 
underscore  the  nature  of  the  approximations  being  made. 
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The  fields  of  the  structure  are  represented  in  the  form  of  the  product  of  a  slowly  varying,  complex 
amplitude  5A„,  periodic  eigenfunctions  e„(x)  and  b„{x),  and  an  exponential  phase  factor. 


Erf(x,t)  SA^iz)  e^{x)  exp  [ii  f  dz'  -  (O^t)]  +  c.c.  (la) 

"■  Jo 

and 

CO 

Brf(x,t)  =  X  j  M„(z)  b^{x)  exp  [/(  k^^iz')  dz'  -  «„0]  +  c.c.  .  (lb) 

n  Jo 


Here  it  is  assumed  that  e„(x)  and  b„(x)  are  the  solutions  of  Maxwell’s  equations  for  the  empty 
structure  corresponding  to  angular  frequency  co„  and  real  axial  wave  number  k^.  In  this  regard,  e„(x) 
and  b„(x)  will  be  periodic  in  axial  distance  with  a  period  equal  to  that  of  the  structure.  To  allow  for 
the  slow  axial  variation  of  the  structure’s  parameters,  the  wave  number  corresponding  to  frequency 
co„  as  well  as  the  eigensolutions  e„(x)  and  b„(x)  will  vary  with  axial  distance.  It  is  assumed  that  such 
variations  are  slow  and  that  it  is  permissible  to  think  of  the  functions  e„(x)  and  b„{x)  as  being  locally 
periodic.  As  Maxwell’s  equations  for  the  cold  structure  are  linear,  we  can  take  the  eigensolutions 
e„(x)  and  b„(x)  to  be  dimensionless.  In  this  case,  the  slowly  varying  amplitude  SA„  has  the  dimensions 
of  a  vector  potential. 


The  sum  over  the  subscript  n  in  Eqs.  (la)  and  (lb)  represents  a  sum  over  frequencies  (0„.  Note  that 
all  the  time  dependence  of  the  fields  is  contained  in  the  exponential  factors;  that  is,  the  amplitudes  5A„ 
depend  only  on  axial  distance.  This  representation  of  the  field  corresponds  to  the  case  in  which  the 
device  is  excited  by  a  signal  consisting  of  a  discrete  set  of  frequencies  belonging  to  the  set  co„.  In 
CHRISTINE,  all  frequencies  are  assumed  to  be  integer  multiples  of  some  minimum  frequency  Aco.  In 
this  case,  the  nonlinear  beating  of  any  two  signals  in  the  device  will  produce  a  signal  at  a  frequency 
that  is  a  member  of  the  set  (0„.  This  beating  could  be  in  the  form  of  self-beating,  which  would 
generate  harmonic  frequencies,  or  it  could  be  in  the  form  of  beating  of  two  signals  of  different 
frequency,  which  would  generate  intermodulation  distortion.  The  restriction  that  members  of  the  set 
co„  are  integer  multiples  of  some  minimum  frequency  Aco  can  be  stated  in  another  way.  The  input 
signal  and  all  other  quantities  of  interest  are  assumed  to  be  periodic  in  time  with  period  T  =  Itt/Ao). 
Section  2.2  discusses  the  implications  of  this  for  the  injection  of  particles  in  the  simulation. 


The  time-averaged  electromagnetic  power  flux  along  the  structure  implied  by  Eqs.  (la)  and  (lb) 
can  be  evaluated  using  Poynting’s  theorem. 


P 


c 

Ik 


2  ^  M„(z) 


where  A,g„  is  the  effective  cross-sectional  area  defined  by  the  relation 


(2) 


^ejfA 


Z  -ielxb^ 


+ 


(3) 


and  z  is  a  unit  vector  in  the  axial  direction.  For  a  structure  with  slowly  varying  parameters,  the 
contribution  to  the  power  flux  from  each  frequency  (i.e.,  each  term  in  the  sum  in  Eq.  (2))  will  be 
independent  of  axial  distance  in  the  absence  of  a  beam  or  attenuation.  Basically,  this  assumes  that  the 
parameters  are  tapered  gently  enough  so  that  the  amount  of  power  reflected  from  a  forward 
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propagating  wave  is  negligibly  small.  With  tapering,  the  quantity  will  vary  with  axial  distance  as 
the  parameters  of  the  structure  vary.  Consequently,  according  to  Eq.  (2),  in  a  structure  with  slowly 
varying  parameters,  the  amplitude  5A„  will  vary  with  axial  distance  even  if  the  power  flux  is  constant. 
To  account  for  this,  we  introduce  a  normalized  field  amplitude 


1/2 


(4) 


where  q  and  m  are  the  charge  and  mass  of  an  electron,  respectively.  The  code  determines  the 
dimensionless  quantity  a„.  This  choice  of  normalization  gives  the  following  expression  for  the  power 
flux: 


P  =  Pflux,2  X I  «n(z)  1^ . 

where 


(5) 


=  1'38S2  X  10’  (6) 

is  a  constant.  A  consequence  of  normalization  is  that  under  the  stated  assumptions,  specifically  the 
negligibly  small  reflection  of  forward  propagating  power,  a„  varies  because  of  attenuation  and 
coupling  to  the  beam,  not  because  of  the  structure’s  slowly  changing  parameters. 

The  complex  amplitude  is  derived  by  following  these  standard  steps.  Expressions  (la)  and  (lb) 
are  inserted  in  Maxwell’s  equations  along  with  the  beam  current.  Ampere’s  law  is  dotted  with  e*(x) 
and  Faraday’s  law  is  dotted  with  b*(x).  The  two  products  are  combined  in  the  same  way  as  one 
forms  Poynting’s  theorem,  and  the  resulting  expression  is  averaged  over  the  spatial  period  of  the 
structure  and  the  temporal  period  T  of  the  fields.  Weak  attenuation  is  introduced  by  the  boundary 
condition  that  the  tangential  electric  field  of  the  radiation  does  not  quite  vanish  on  the  structure  walls. 
Thus,  the  resulting  equation  for  the  amplitude  of  the  n'*'  spectral  component  is 

^  +  a„(z))  a„  =  y^72  (J  ^^P  t ^zn(z')  dz'  -  0)^t)]j  ,  (7) 

where  an(z)  is  the  attenuation  in  Nepers/cm  and  is  frequency  dependent,  I  a  =  mc^/q,  j  is  the  beam 
current  density,  and  the  angular  brackets  denote  averages  over  the  temporal  period  of  the  radiation 
and  the  spatial  period  of  the  structure.  The  right-hand  side  of  Eq.  (7)  is  further  refined  in  Section  2.2 
after  a  discussion  of  the  equations  of  motion. 
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2.2  Equations  of  Motion 

In  CHRISTINE,  beam  electrons  are  treated  as  annular  discs  of  charge  of  outer  radius  and  inner 
radius  which  are  constrained  to  move  in  only  the  axial  direction.  One  expects  this  approximation 
to  be  appropriate  when  the  electrons  are  confined  by  either  a  strong  axial  magnetic  field  or  a  system 
of  periodic  focusing  magnets.  The  motion  of  the  electrons  will  be  one  dimensional  if  the  frequency 
of  the  amplified  signal  is  much  smaller  than  the  frequency  of  particle  oscillations  in  the  relevant 
confining  field  (the  gyrofrequency  or  the  betatron  frequency).  The  disc  approximation  will  be 
appropriate  if  the  axial  electric  field  of  the  signal  does  not  vary  appreciably  over  the  radial  extent  of 
the  electron  beam. 

The  rate  of  change  of  electron  energy  is  expressed  in  terms  of  the  rate  of  change  of  the  relativistic 
factor  7  =  1/  >J\  —v^tc^  , 


dr  ^  q 

dz  mc^ 


{  ^  ■  {^rf  +  ^sc  +  beam  ' 


(8) 


Here,  Erf  is  the  electric  field  of  the  structure  field  given  by  Eq.  (la),  Esc  is  the  AC  space  charge  field 
discussed  in  Section  2.3,  and  E^c  is  a  steady-state  axial  electric  field,  if  present,  that  changes  the 
energy  of  beam  electrons  as  they  travel  through  the  structure.  A  steady-state  axial  field  will  be  present 
if  the  DC  space  charge  depression  varies  with  axial  distance  (e.g.,  if  the  parameters  of  the  structure  are 
tapered).  Equation  (8)  describes  the  Lagrangian  rate  of  energy  change  of  a  particle  as  it  travels  down 
the  interaction  region.  The  independent  coordinate  is  the  axial  location  of  an  electron.  To  determine 
the  arrival  time  t{z)  of  an  electron  at  a  particular  axial  location,  we  solve 


dt  _  1 

dz  ~  vlf)  ’ 


(9) 


where  the  axial  velocity  of  an  electron  is  given  in  terms  of  the  relativistic  factor  y  and  the  pitch  factor 
0is 


yl  j 


(10) 


The  pitch  factor  0  is  the  ratio  of  the  transverse  component  to  the  total  momentum  on  injection  and  yo 
is  the  electron’s  relativistic  factor  on  injection.  The  model  assumes  that  0  «1,  i.e.,  small  enough  so 
that  the  motion  is  one  dimensional  but  large  enough  that  there  can  be  a  significant  spread  in  axial 
velocities.  In  CHRISTINE,  it  is  assumed  that  the  pitch  factor  has  a  Gaussian  distribution  and  that  the 
RMS  width  of  the  distribution  can  be  specified  (see  Section  3).  The  initial  conditions  for  the 
equations  of  motion  are  the  following:  particles  are  injected  with  a  specified  relativistic  factor  pitch 
factor  0,  and  entrance  time  t(0).  The  number  and  distribution  of  injected  particle  entrance  times  are 
discussed  later  in  this  section. 
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CHRISTINE  does  not  solve  the  quantity  t(z)  directly  using  Eq.  (9).  Instead,  it  determines  the  phase 
relative  to  the  n'’'  signal.  This  phase  is  defined  as 


¥n  = 


(11) 


where  v^o  is  the  initial  axial  velocity  of  a  particle  with  pitch  factor  zero.  It  is  evaluated  taking  voltage 
depression  into  account.  The  evolution  of  this  phase  is  then  given  by 


dWn 

dz 


ty„(l/v^-l/v,). 


(12) 


Introducing  the  phase  Xj/n  into  the  expression  for  the  radiation  field  (Eq.  (la))  gives 


dy 

dz. 


rf 


Re  |2i  2  a„(z)  e2(n,z) 


(13) 


for  the  contribution  to  the  rate  of  change  of  ydue  to  the  structure  fields.  In  Eq.  (13),  the  quantity 
e2(n,z)  is  defined  as 


e2(.n,z) 


{z  •  en)be 

.Ml 


exp 


JO 


(z')  -  ®„/Vj^)r/z'] 


(14) 


where  the  angular  bracket  denotes  average  over  the  radial  extent  of  the  beam  and  over  one  axial 
period  of  the  structure  fields. 


The  above  quantity  also  appears  effectively  in  Eq.  (7)  for  the  complex  field  amplitude.  To  see  this, 
we  first  represent  the  current  density  in  the  form  of  a  sum  over  moving  charge  discs 


5{z.  -Zj{t)) 


0,  r  <ri,i  ) 

1 ,  r^f<r  <  ri,A, 
0,  r>  j 


(15) 


where  the  sum  is  over  all  the  electrons.  Equation  (7)  calls  for  an  average  over  one  spatial  period  of 
the  structure  and  one  temporal  period  of  the  fields.  Due  to  the  time  periodicity  of  all  quantities,  this 
average  may  be  evaluated  by  integration  over  the  shaded  region  shown  in  Fig.  1 . 
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Fig.  1  —  Characteristics  in  the  r  vs  z  plane  for  electron  trajectories.  The  shaded  region 
corresponds  to  one  spatial  period  of  the  structure  and  one  temporal  period  of  the  fields. 


To  evaluate  the  average,  the  shaded  region,  a  parallelogram,  is  preferable  to  a  rectangle,  with  upper 
and  lower  boundaries  at  fixed  t  because  it  contains  the  orbits  of  a  group  of  electrons  that  all  pass 
through  both  lateral  boundaries.  The  average  in  Eq.  (7)  may  be  evaluated  by  first  integrating  over 
time: 

f  v,,/0  S{z  -  zjit))  exp  [  -  i{  f'  k^iz')  dz'  -  co^t)-] 

t  ^  Jo 


. 


dz  e*(jc)  •  z 


exp  [  -  /(J^  dz'  - 


(16) 


where  the  sum  over  particles  represents  a  sum  over  the  group  of  particles  entering  the  interaction 
region  during  one  period  of  time  T  as  depicted  in  Fig.  1  and  tjiz)  is  the  time  of  flight  for  particle  j'io 
the  point  z.  The  number  of  particles  that  will  contribute  to  the  sum  is  the  number  that  entered  during 
a  time  T,  that  is,  Tl/q,  where  7  is  the  beam  current.  Now  the  integral  over  z  can  be  carried  out.  Because 
the  exponent  is  evaluated  along  a  particle  trajectory,  it  can  be  regarded  as  slowly  varying  over  a 
distance  given  by  the  structure  period  A//.  Thus,  the  average  over  z  reduces  to  an  average  of  the  axial 
component  of  the  structure  field.  Including  the  radial  dependence  of  the  current  density  implied  by 
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the  disc  model,  we  see  that  Eq.  (7)  calls  for  precisely  the  same  average  of  the  structure  fields  as  does 
the  equation  of  motion.  The  result  is  that  the  average  in  Eq.  (7)  produces 

(ji  ^ 

y4 

where  7  is  the  total  beam  current  and  the  angular  brackets  now  signify  average  over  particles  entering 
the  interaction  region  during  the  time  interval  T.  Since  the  beam  current  enters  in  the  form  of  a  ratio 
I/I  A,  the  beam  current  and  Ia  can  be  evaluated  in  any  system  of  units.  In  CHRISTINE,  I  is  evaluated 
in  amperes  and  Ia=  1.7  x  10"  amps.  It  is  clear  from  Eqs.  (13)  and  (15)  that  the  relevant  information 
concerning  the  structure  fields  is  contained  in  the  complex  coupling  function  e2(n,z).  The  amplitude 
of  this  quantity  is  related  to  the  axial  impedance.  In  particular, 

\p  -  1^^  '^n)beam\  _  ^ 


where  K  is  the  coupling  impedance  in  ohms.  The  phase  of  the  coupling  function  (Eq.  (14))  is 
determined  by  the  phase  velocity  of  the  structure  fields. 

The  initial  conditions  on  the  phase  are  determined  by  the  distribution  of  particle  entrance 
times.  Basically,  entrance  times  must  be  prescribed  over  the  temporal  period  of  the  fields  T.  If  the 
electron  beam  is  unmodulated,  the  entrance  times  are  uniformly  distributed  over  this  time  interval.  In 
the  case  of  a  premodulated  beam  with  ballistic  bunching,  the  initial  entrance  times  are  distributed 
according  to  the  formula 


7(Z  =  0)  =t' 


'L^cos(m,t'  +0BJ 

n  =  1 


(19) 


where  the  times  f^are  uniformly  distributed  over  the  interval  T.  In  Eq.  (19),  QB„  and  0Bfi  are  the 
amplitude  and  phase  of  the  modulation  signal  for  frequency  (On-  In  the  case  of  premodulation  at  a 
single  frequency  oo^,  the  ratio  of  average-to-peak  beam  current  implied  by  Eq.  (19)  is  1  \ . 

Gated  emission  modulation  presents  two  options.  Either  it  is  assumed  that  the  current  is  on  or  off, 
giving  a  “top  hat”  dependence  to  the  current  on  time,  or  it  is  assumed  that  the  current  has  a  “sin^” 
dependence  on  time.  In  these  cases,  the  entrance  times  are  uniformly  distributed  over  a  subinterval  of 
the  period  of  the  bunching  signal.  The  duration  of  the  subinterval  is  specified  by  giving  the  ratio  of 
the  average  current  to  the  peak  current.  In  the  “sin^”  case,  each  entering  particle  is  given  a  weight 
between  0  and  2  according  to  the  “sin^”  distribution. 

2.3  The  Space  Charge  Field 

This  section  describes  the  procedure  for  modeling  the  space  charge  electric  field  shown  in  Eq.  (8). 
The  effect  of  this  field  is  to  resist  the  formation  of  electron  bunches  and  consequently  to  give  rise  to 
collective  oscillations  of  the  beam  at  a  modified  plasma  frequency.  The  natural  frequency  of  these 
oscillations  depends  on  the  beam  density,  its  radial  profile,  and  the  configuration  of  the  metallic 
structure  surrounding  the  beam.  We  assume  that  the  basic  form  of  the  field  can  be  written  as 
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^  — rr . 

'■wj 


(20) 


where  /?„'  is  an  unknown  coefficient  that  will  be  determined  subsequently.  Were  this  the  only  field 
acting  on  the  particles,  the  linearized  bunching  factor  would  oscillate  in  time  with  an  angular 
frequency  where 


,2  _  4g//? 


This  result  is  obtained  by  taking  the  small  signal  limit  of  Eqs.  (8),  (10),  and  (12)  and  assuming  for 
simplicity  that  the  beam  is  cold  (0  =  0  in  Eq.  (10)).  The  factor  R'  describes  the  reduction  in  the 
collective  oscillation  frequency  accounting  for  the  distribution  of  fields  outside  the  beam. 

CHRISTINE  uses  the  sheath  helix  model  to  calculate  the  effects  of  the  structure  in  the  space 
charge  field.  This  approach  differs  from  others  in  that  no  attempt  is  made  to  separate  fields  into  their 
electromagnetic  and  electrostatic  components.  This  also  allows  unambiguous  definition  of  the  correct 
coefficient  for  the  space  charge  field  to  include  in  the  equation  of  motion.  The  space  charge  field 
that  we  seek  has  a  time  dependence  similar  to  that  of  the  structure  fields  (i.e.,  it  is  periodic  in  time 
with  period  T  and  thus  includes  frequencies  belonging  to  the  set  cOn).  In  general,  more  frequencies 
are  needed  to  correctly  model  the  space  charge  term  than  to  model  the  structure  fields.  This  is 
because,  in  the  nonlinear  regime,  the  bunched  current  density  of  the  beam  has  a  high  harmonic 
content.  Frequency  co„  is  generally  accompanied  by  spatial  wave  number  These  spectral 

components  of  the  current  tend  to  excite  modes  of  the  structure  over  that  range  of  frequencies  where 
the  dispersion  relation  of  the  structure  is  linear.  However,  as  the  excitation  of  the  space  charge  field 
does  not  rely  on  resonance  with  the  structure  fields,  the  harmonic  content  of  the  current  density 
produces  appreciable  space  charge  fields  over  a  greater  range  of  frequencies. 

The  approach  taken  by  CHRISTINE  is  to  assume  that  the  beam  current  density  is  known  and  to 
then  calculate  the  total  (‘structure’  plus  ‘space  charge’)  axial  electric  field  at  the  beam  that  appears 
in  response  to  the  given  current  density.  From  the  total  field,  one  can  then  extract  the  appropriate 
factor  R'.  We  assume  that  the  form  of  the  current  density  corresponds  to  an  ensemble  of  moving 
annular  discs  oscillating  in  time  with  a  frequency  %  and  an  as  yet  unspecified  wave  number  k^.  In 
order  to  connect  the  present  calculation  with  the  model  of  the  current  used  in  Eq.  (7),  we  write 


r  <rt,i  ) 

,  ryi<r  <  ri,A+c.c.  , 
r>  J 


(22) 


where  and  are  the  outer  and  inner  radii  of  the  beam  discs,  and  the  angular  brackets  represent 
the  ensemble  average  over  particles.  Inserting  Eq.  (22)  and  its  corresponding  charge  density 
perturbation  determined  by  continuity  into  Maxwell’s  equations  results  in  the  following  radial 
differential  equation  for  the  complex  amplitude  of  the  axial  electric  field  with  frequency  cOn  and 
wavenumber 
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4/ 


0 ,  r  <  I 

1  ,  ri,i<r  <  ri,A  , 

r>  ri,„  I 


(23) 


where 


(24) 


We  solve  Eq.  (23)  for  r  <  r„,  the  radius  of  the  sheath  helix.  In  addition,  we  solve  the  similar  (but 
homogeneous)  equation  for  the  axial  magnetic  field  of  the  TE  component  of  the  radiation  for  r  <  r„. 
Outside  the  helix  (i.e.,  r  >  r„),  we  solve  similar  homogeneous  equations  for  the  TE  and  TM 
components  of  the  field,  but  assume  that  a  dielectric  with  relative  permitivity  e  is  present.  The 
following  boundary  conditions  are  imposed:  EQ(r)  vanishes  at  r  =  rjv  >  r„;  E^(r)  vanishes  at  r=  rv>  r„ 
(where  >  Vy  ),  simulating  the  effect  of  vanes;  and  the  components  of  the  electric  field  parallel  to 
the  helix  direction  and  the  surface  current  perpendicular  to  the  helix  direction  vanish  at  the  helix. 
Finally,  we  average  the  axial  electric  field  over  the  annular  discs.  The  resulting  expression  for  this 
averaged  field  is 


^  beam 


2i(0„ 


c^D{k^,0)„) 


Ail 


(25) 


The  quantity  D(kz,(On)  in  the  denominator  of  the  first  term  in  Eq.  (25)  is  the  vacuum  dispersion 
function  for  the  sheath  helix: 


D(A:,,a>J 


c2  [KTf^lQiKrf,) 


K^h 


'^H^h 


KThloiKT^) 


(26) 


where  k^  is  the  helix  wave  number,  and  Cy  are  linear  combinations  of  Bessel  functions 


and 


CyiN^h)  =Io(^^h)KQ(N^^)-KQ(K^h)lQ(K^^)  , 


(27a) 

(27b) 


(27c) 
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In  Eq.  (25),  quantities  H  and  R  are  defined  as  follows: 


and 


H  = 


R  =  l  +  - 


2  2 
rbo-fti 


r/iM;: 


rbo  Ko  (»6«  )  -  1  ('"■)] 


(28a) 


(28b) 


The  first  term  in  Eq.  (25)  is  essentially  the  contribution  of  the  structure  field  to  the  axial  electric  field 
at  the  beam,  and  the  second  term,  involving  R,  is  essentially  the  space  charge  field.  The  quantity  R  is 
the  space  charge  reduction  factor,  and  it  is  an  interesting  exercise  in  Bessel  function  identities  to  show 
that  R  vanishes  as  rt,i  r^o- 

To  determine  the  space  charge  field  to  be  inserted  in  Eq.  (8),  we  first  write  it  in  the  form 


{^■^sc)bea.=-l.^  (^2 

n  ^n\'bo  ^bi) 


(29) 


with  R'n  to  be  determined.  Assuming  all  fields  vary  as  exp(ikzz),  the  total  field  implied  by  Eqs.  (29) 
and  (7)  is 


{z-E^  +  z-E,,) 


2m 


sc/  beam  Y  c[k^  - 


I  (e-'>) 


2 

z-en 

^eff,n 


Ail 


+  C.C.. 


(30) 


As  can  be  seen,  Eqs.  (25)  and  (30)  have  similar  behavior  with  respect  to  their  dependence  on  k^.  They 
both  diverge  as  k^  k^n,  where  k^n  is  the  solution  of  the  vacuum  dispersion  relation  at  frequency  (On- 
We  pick  the  space  charge  factor  Rn'lo  give  agreement  between  Eqs.  (25)  and  (30)  when  both  are 
expanded  about  k^  =  k^n-  To  lowest  order,  the  coefficient  of  the  singular  term  determines  the 
impedance  in  the  sheath  helix  model 


^  ^  (On 

377  Q  c  dD{kp(0^  Idk^ 


(31) 


In  next  order,  matching  the  coefficient  of  the  constant  term  gives 
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where 


and 


R'  =  R- 


4c^{D:f 


d:  = 


d{D/H^) 


dk. 


dk. 


(32) 


(33a) 


(33b) 


The  difference  between  R  and  ^'can  be  interpreted  as  the  effect  of  the  sheath  helix  on  the  space 
charge  reduction  factor.  In  particular,  if  the  helix  is  replaced  by  a  conducting  tube,  no  structure  field 
can  resonate  with  the  beam,  and  the  space  charge  reduction  factor  is  simply  R.  With  the  helix  present, 
some  of  the  space  charge  field  is  able  to  penetrate  outside  the  helix,  and  the  reduction  factor  is 
modified. 


Thus,  the  final  equation  of  motion  is 


dz 


8/i 

h{fbo-rbi 


c_K 


+ 


beam  ’ 


(34) 


where  the  separate  sums  over  n  and  n'are  over  frequencies  belonging  to  the  set  cOn-  In  CHRISTESE, 
the  number  of  frequencies  kept  to  evaluate  the  space  charge  term  is  independent  of  the  number  kept 
to  evaluate  the  structure  fields.  In  particular,  a  range  of  frequencies  is  specified  in  the  input  list  for  the 
radiation  field  as  is  the  number  of  space  charge  harmonics.  In  computing  the  space  charge  factor  Rn' 
for  each  frequency  harmonic,  the  code  first  checks  to  see  if  the  frequency  is  in  the  range 
corresponding  to  the  structure  fields.  If  it  is,  the  reduction  factor  is  calculated  according  to  Eq.  (32). 
If  the  space  charge  harmonic  frequency  is  out  of  the  range  specified  for  the  structure  fields,  the 
reduction  factor  is  taken  to  be  just  R  and  calculated  for  a  spatial  wave  number  k^  =  (On/v^o- 

3.  NUMERICAL  IMPLEMENTATION 

This  section  outlines  the  numerical  algorithm  used  to  solve  the  basic  equations.  At  various  points 
in  the  discussion,  physical  variables  from  Section  2.2  are  identified  along  with  their  FORTRAN 
equivalents  in  the  program,  which  appear  in  brackets. 

The  main  numerical  calculation  in  the  program  is  the  integration  of  Eqs.  (12),  (15),  and  (34).  This 
is  done  by  a  fourth  order  Runga-Kutta  routine  in  subroutine  ‘rk’.  This  subroutine  calls  subroutine 
‘eqm’  to  evaluate  the  derivatives  of  each  function.  The  following  variables  are  integrated:  the 
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complex  field  amplitude,  an(z),  [as(nn,j)];  the  relativistic  factor,  y(z),  [gamma(kk)];  and  the  sine  and 
cosine  of  the  particle  phase. 

Because  the  complex  exponential  of  the  phase  of  each  particle  with  respect  to  perturbations  at 
each  frequency  is  ultimately  needed,  i.e.,  Z„(zj  =  exp(  -iy/niz)),  [zc(nn,kk)],  the  following  procedure  is 
used.  The  phase  is  defined  in  Eq.  (11),  and  its  complex  exponential  is  used  in  Eqs.  (13),  (15),  and 
(34).  The  values  of  the  Zn  are  not  independent  because  the  frequencies  are  equally  spaced. 
CHRISTINE  calculates  the  sine  [s(kk)]  and  cosine  [c(kk)]  of  the  phase  for  mode  0,  which 
corresponds  to  the  central  frequency  (Oq: 


=  0}q{1  /  v^-  1  /  V,)  cos  y/Q , 


dz 

d  cos  y/Q 
dz 


=  -(Oq(\ /v^-l /v^)sini(fQ, 


(35a) 


(35b) 


and  the  sine  [sd(kk)]  and  cosine  [cd(kk)]  of  the  phase  corresponding  to  the  difference  frequency 
Aco: 

d  sin  Aiff 


dz 


=  Ao)  (1  /  - 1  /  v^)  cos  Ay/, 


d  cos  Aw  A  ^  I  \  ■  A 

- ^  =  -Aci){\lv^-\lv^smAy/. 

The  following  quantity  is  then  formed  in  subroutine  ‘eqm,’ 

Z^-e  "t'"  =  (cos  y/Q  +  i  sin  y/Q)(cos  Ay/  +  i  sin  Ay/Y 


(36a) 


(36b) 


(37) 


In  addition,  the  above  quantity  is  raised  to  the  m”'  power  to  calculate  the  m'*  harmonic  of  the  space 
charge  field  for  each  frequency  [zcm(nn,kk)]. 


=  (ZJ" . 


(38) 


Some  input  parameters  may  overlap  because  the  m**  harmonic  of  the  space  charge  field  for  one 
frequency,  nl,  coincides  with  a  higher  frequency  n2.  This  will  occur  \i  m  {(Oq  A  nl  Am)  =  (Oq  +  nl 
Aco.  Such  possibilities  are  detected  in  subroutine  ‘initparam’,  where  a  table  [iten(nn,nh)]  is  generated. 
If  the  element  of  the  table  for  a  particular  n  and  m  value  is  unity,  then  the  corresponding  Z^^rn  is  kept 
in  the  space  charge  term.  If  the  value  is  zero,  the  term  is  not  kept. 
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3.1  Subroutines 

The  following  paragraphs  list  and  describe  the  functions  of  the  various  subroutines  used  in 
CHRISTINE. 

3.1.1  Subroutine  ‘initparam’ 

This  subroutine  performs  a  number  of  preliminary  calculations  prior  to  integrating  the  wave 
equation  and  equations  of  motion.  Among  these  are 

1)  introducing  various  physical  constants, 

2)  determining  the  frequencies  present  in  the  RF  spectrum, 

3)  determining  the  total  number  of  particles, 

4)  invoking  either  subroutine  ‘S_HELIX’  to  calculate  the  phase  velocity  and  impedance  for 
the  sheath  helix  model,  or  subroutine  ‘READ_DAT’  to  read  in  this  information, 

5)  generating  the  axial  grid  and  the  interpolation  factors  for  the  profile  of  tapered  quantities, 

6)  calculating  the  accelerating  gradient  [dvbeam(i)], 

7)  calculating  the  attenuation  profile, 

8)  generating  the  distribution  of  electron  velocities  for  a  warm  beam, 

9)  calculating  the  space  charge  coefficient  for  each  frequency,  and 

10)  initializing  the  amplitude  of  the  injected  signal  at  each  frequency. 

3.1.2  Subroutine  ‘expon’ 

This  subroutine  calculates  the  coupling  factor  e2(n,z)  [e2(nn,j)],  which  is  defined  in  Eq.  (14).  It 
requires  the  phase  velocity  and  coupling  impedance  for  each  frequency. 

3.1.3  Subroutine  ‘phases’ 

This  subroutine  initializes  the  phases  of  the  entering  simulation  particles  as  detailed  in  Section  2. 
For  a  gated  emission,  weights  are  also  given  to  the  entering  particles  to  describe  the  time  dependence 
of  the  beam  current. 

3.1.4  Subroutine  ‘rk’ 

This  subroutine,  the  Runga-Kutta  integrator,  integrates  the  equations  of  motion  and  the  wave 
equation.  It  invokes  subroutine  ‘eqm’  to  calculate  the  rate  of  change  of  each  quantity  to  be  varied.  In 
addition,  this  subroutine  tabulates  and  writes  the  particle  data  to  files.  This  information  is  discussed 
more  completely  in  Section  4.1. 

3.1.5  Subroutine  ‘srs’ 

This  subroutine  invokes  subroutine  ‘rk’  and  it  tracks  the  power  balance. 
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3.1.6  Subroutine  ‘eqm’ 

This  subroutine  computes  the  right-hand  side  of  Eqs.  (15),  (34),  (35),  (37),  and  (38).  It  is  invoked 
by  subroutine  ‘rk’. 

3.1.7  Subroutine  ‘S_HEL1X’ 

This  subroutine  solves  the  dispersion  relation  for  the  sheath  helix  model  by  using  Newton’s 
method.  It  also  determines  the  coupling  impedance  and  space  charge  correction  factor  caused  by 
fields  outside  the  helix. 

3.1.8  Subroutine  ‘READ_DAT’ 

This  subroutine  reads  in  the  parameters  of  the  structure  as  an  alternative  to  using  the  sheath  helix 
model.  Section  4  describes  the  format  for  inputting  data. 

3.2  Input  Namelists 

Input  for  CHRISTINE  is  provided  through  ‘tiniS.nml’  (Table  1),  a  file  that  contains  a  sample  set 
of  namelists  corresponding  to  a  specific  research  tube.  This  input  file  is  used  as  a  basis  for  the  set  of 
example  runs  presented  in  Section  4.  The  various  namelists  are  summarized  in  the  following 
paragraphs. 


Table  1  —  CHRISTINE  Input  File  ‘tiniS.nml’ 

curbeam  -  beam  current  (in  Amps) 
tht_sprd  -  maximum  injection  angle  in  degrees 
controls  effective  axial  velocity  spread 
vbeam  -  beam  voltage  (in  Volts) 
rbo  -  outer  beam  radius  (in  cm) 
rbi  -  inner  beam  radius  (in  cm) 

shotn  -  shot  noise  in  beam  (shotn  =  1 .  is  the  nominal  value) 
dc_sc  -  if  true,  DC  space  charge  depression  is  applied  to  beam  energy 

$beam 

curbeam  =  .17, 
tht_sprd  =  0., 
vbeam  =  2.84e3, 
rbo  =  .036, 
rbi  =  .03, 
shotn  =  0., 
dc-sc  =.false. 

Send 

delfreq  -  resolution  in  frequency 
freqO  -  central  frequency 
freq_max  -  maximum  frequency 

freq_min  -  minimum  frequency 

pow_in  -  array  of  input  powers  [watts] 
phase_in  -  array  of  input  phases  in  radians 
QB  -  array  of  prebunching  factors 
Phase-B  -  array  of  prebunching  phases 

gated_I  -  if  true,  modulation  of  current _ 
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Table  1  —  CHRISTINE  input  file  ‘tiniS.nml’  (continued) 

sin_2  -  if  true,  sin^2  time  dependence,  otherwise  top  hat 
Avel_pkl  -  average  current  divided  by  peak  current 

$rad 

freqO  =  5.e9, 
delfreq  =  5.e9, 
freq_max  =  5.e9, 
freq_min  =  5.e9, 
pow_in  =  2.4e-03, 
phase_in  =  0., 

QB  =.0, 
phase_B  =  0., 
gated_I  =  .false., 
sin_2  =  .true., 

Avel_pkl  =  .60, 

Send 

nz  -  number  of  axial  points 

naO  -  number  of  velocity  angles 

fntO  -  number  of  entrance  times  per  mode 

npO  -  number  of  phases 

nsc  -  number  of  space  charge  harmonics 

Snum 

nz  =  300, 
naO  =  1, 
fntO  =1.5, 
npO  =  137, 
nsc  =  1, 

Send 

zint  -  interaction  length 

zr(i)  -  axial  points  where  parameters  are  specified. 

All  quantities  are  piecewise  continuous  between 
these  points,  the  first  point  should  be  zero  and 
the  last  should  be  greater  than  or  equal  to  zint. 
dvbeam(i)  -  acceleration  potential 
use_sh  -  if  true,  uses  sheath  helix  model 

Sstruct 

zint  =  9.5758, 
zr  =  0.,  10., 
dvbeam  =  2*0. 
use_sh  =  .true.. 

Send 

RH(i)  -  array  of  helix  radii 

H_LMDA(i)  -  array  of  helix  periods 

RW(i)  -  array  of  wall  radii 

EPS_R(i)  -  array  of  dielectric  constants 

RV(i)  -  array  of  vane  radii  (if  no  vanes,  set  this  equal  to  RW) 

use_cor  -  if  true,  adds  correction  to  space  charge  term 
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Table  1  —  CHRISTINE  input  file  ‘tiniS-nml’  (continued) 

$helix_par 
RH  =  2*.  12446, 

H_LMDA  =  2*.080137, 

RW  =  2*.2794, 

EPS_R  =2*1.75, 

RV=  2*  .2794, 
use_cor  =  .false. 

Send 

za(i)  -  array  of  z  points  for  specification  of  attenuation 
att(i)  -  array  of  attenuation 

dB_cm  -  if  true,  att  in  dB/cm,  otherwise  dB/wavelength 
f_scale  -  attenuation  proportional  to  (f/freqO)**f_scale 

$loss_par 
za  =  0.,10., 
att  =  2*0., 
dB_cm  =  .true., 
f_scale  =  0. 

Send 

nt_plot  -  number  of  times  for  beta  vs  z  plots 
ps_z  -  array  of  axial  positions  for  phase  space  plots 

Sps_plot 
nt_plot  =  3, 
ps_z  =  5.,7.,9., 

Send 

To  read  in  structure  properties 
npts  -  number  of  data  points 
f_in(i)  -  array  of  frequencies 
bp_in(i)  -  array  of  phase  velocities 
z_in(i)  -  array  of  impedances 

RH_eff  -  effective  helix  radius  for  calculating  space  charge 
fract_bp  -  modify  above  phase  velocities  at  points  zr(i) 
fract_Z  -  modify  above  impedances  at  points  zr(i) 
use_dat  -  if  true,  read  in  parameters  from  my  .bp  and  my.zc 

Sdis_dat 
npts  =  15, 
f_in  = 

I.e09,2.e09,3.e09,4.e09,5.e09,6.e09,7.e09,8.e09,9.e09,10.e09,ll.e09,12.e09,13.e09,14e09,15.e09, 
bp_in  = 

0. 1 0790,0. 1 0530,0. 10160,0.097790,0.094480,0.092000,0.090270,0.089090,0.088300,0.087770,0.0 
87400,0.087140,0.086960,0.086840,0.086750 
z_in  = 

236.40,212.30,172.00,123.20,78.920,47.030,27.090,15.430,8.8040,5.0490,2.9170,1.6970,0.99390,0 

.58570,0.34710 

fract_bp  =2*1.,.89,.79 
fract_Z  =  2*1.,.84,.68 
RH_eff  =  .12446, 


CHRISTINE 


17 


Table  1  —  CHRISTINE  input  file  ‘tiniS.nml’  (continued) 

use_data  =  .true., 

Send 

scan_p  -  parameter  scanning  activated  if  .true. 
n_hann  -  number  of  harmonics  in  radiation  field 
(  freq_max  =  (n_harm+l)*freqO  ) 

varX  -  variable  to  be  scanned 
varY  -  variable  to  be  scanned 
possibilities  are: 

curbeam 

vbeam 

tht_sprd 

rbo 

rbi 

ffeqO 

pow_in  (1) 
pow_in  (2) 
phase_in  (1) 
phase_in  (2) 

QB 

Avel_pkl 

nx  -  number  of  values  of  varX 

ny  -  number  of  values  of  varY 

varX_min  -  minimum  value  of  varX 
varX_max  -  maximum  value  of  varX 
varY_min  -  minimum  value  of  varY 
varY_max  -  maximum  value  of  varY 
log_X  -  if  true,  varX  scanned  logarithmically 
log_Y  -  if  true,  varY  scanned  logarithmically 

Sscan 

scan_p  =  .true. 
n_hann  =  0, 
varX  =  ‘vbeam’, 
nx  =  25, 

varX_min  =  2.0e3, 
varX_max  =  3.5e3, 
log_X  =.false. 
varY  =  ‘freqO’, 
ny  =1, 

varY_min  =  5.e9, 
varY_max  =  5.e9, 

Iog_Y  =.false. 

Send 
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3.2.1  Namelist  -  beam 

These  are  the  physical  parameters  of  the  beam.  For  a  cold  beam,  ‘tht_sprd’  =  0.  To  simulate  a 
warm  beam,  set  ‘tht_sprd’  to  the  RMS  width  A6  (in  degrees)  of  the  distribution  of  injected  angles  of 
the  beam  electrons.  CHRISTINE  assumes  the  distribution  is  Gaussian: 


P(0) 


26 
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(39) 


The  number  of  injection  angles  used  is  controlled  by  the  variable  ‘naO’  appearing  in  namelist 
‘num’.  Injection  angles  are  chosen  so  that  the  particles  have  equal  weights  (i.e.,  the  values  of  the 
cumulative  distribution  function  are  equally  spaced): 


naO 


naO  +1/2 


-T j,  i  =\,  naO  . 


(40) 


The  beam  is  assumed  to  be  annular  with  inner  and  outer  radii  ‘rbi’  and  ‘rbo’  such  that  ‘rbi’  < 
0.95  ‘rbo’.  Otherwise,  the  space  charge  term  is  not  calculated  accurately.  For  a  solid  beam,  set  ‘rbi’ 
to  0. 


The  variable  ‘shotn’  controls  the  noise  on  the  injected  beam.  Basically,  each  injected  particle  is 
given  a  small  random  component  to  its  weight.  With  ‘shotn’  =  0,  the  random  component  is  zero  and 
the  simulation  is  noiseless.  With  ‘shotn’  =  1,  the  random  component  is  calculated  to  correspond  to  the 
shot  noise  produced  by  uncorrelated  electrons  at  the  specified  beam  current.  This  is  probably  an 
overestimate  of  the  noise  level  produced  by  space  charge  limited  guns.  The  correct  value  of  noise 
depends  on  features  of  the  gun  not  modeled  in  CHRISTINE. 

The  beam  energy  will  be  depressed  by  the  DC  space  charge  potential  when  the  logical  variable 
‘dc_sc’  is  set  to  .true..  The  amount  of  depression  is  given  by  the  formula 
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where  r„  is  the  helix  radius,  which  is  assumed  to  be  at  ground  potential. 
3.2.2  Namelist  -  rad 


This  namelist  specifies  aspects  of  the  radiation.  The  first  four  variables  determine  the  properties  of 
the  spectrum  of  the  RF  fields.  The  spectrum  is  assumed  to  consist  of  frequencies  spaced  equally 
about  ‘freqO’  by  an  amount  ‘delfreq’.  The  quantity  ‘freqO’  should  be  an  integer  multiple  of 
‘delfreq’.  If  it  is  not,  the  code  adjusts  the  value  of  ‘delfreq’  so  that  it  is.  The  spectrum  then  consists 
of  all  frequencies  between  ‘freq_min’  and  ‘freq_max’  that  differ  from  ‘freqO’  by  an  integer 
multiple  (positive  or  negative)  of  ‘delfreq’.  Arrays  ‘powjn’,  ‘phase_in’,  ‘QB’,  and  ‘phase_B’ 
control  the  signal  injected  into  the  interaction  regions.  The  elements  of  the  array  correspond  to  the 
frequencies  in  the  spectrum  going  from  smallest  to  largest.  Thus,  some  care  needs  to  be  exercised  in 
assigning  numerical  values  to  these  variables.  One  has  to  first  compute  which  frequencies  will  be 
present  in  the  simulation  based  on  the  values  of  ‘freqO’,  ‘delfreq’,  ‘freq_min’,  and  ‘freq_max’,  and 
then  enter  the  numerical  values  for  the  signal  parameters  in  the  correct  order.  The  definition  of  QB 
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and  its  phase  ‘phase_B’  can  be  found  in  Eq.  (19).  The  phases  ‘phasejn’  and  ‘phase_B’  are  defined 
with  respect  to  some  arbitrary  reference  time.  In  principle,  only  the  relative  phases  of  the  input  signals 
affect  the  output  signal.  However,  as  discussed  in  the  next  paragraph,  one  of  the  diagnostics  is  keyed 
to  the  reference  time. 

An  alternate  way  of  exciting  a  signal  in  the  structure  is  to  modulate  the  beam  current.  This  is  done 
by  using  the  logical  input  gated J  =  .true..  With  this  option,  only  one  frequency  and  possibly  its 
harmonics  can  be  used  in  the  simulation.  In  other  words,  ‘delfreq’  is  set  equal  to  ‘freqO’.  The 
temporal  shape  of  the  modulated  current  is  controlled  by  the  logical  variable  ‘sin_2’.  With  sin_2  = 
.true.,  the  current  pulse  has  a  sine-squared  temporal  waveform,  otherwise  it  is  a  top  hat.  The  duration 
of  the  current  pulse  is  controlled  by  the  variable  ‘Avel_pkl’,  which  is  the  ratio  of  the  average  current 
to  the  peak  current  in  the  pulse.  The  variable  ‘curbeam’  always  represents  the  average  current. 

3.2.3  Namelist  -  num 

This  namelist  controls  the  basic  numerical  parameters  in  the  integration  scheme  in  CHRISTINE. 
The  number  of  axial  points  required,  ‘nz’,  depends  on  the  number  of  exponentiations  of  the  signal 
as  it  passes  through  the  interaction  region.  A  value  of  200  is  probably  a  safe  minimum;  if  a  much 
larger  value  is  required,  the  results  should  be  tested  for  significant  change.  The  number  of  injected 
particles  is  controlled  by  the  variables  ‘naO’,  ‘fntO’,  and  ‘npO’.  As  mentioned,  ‘naO’  controls  the 
number  of  injected  pitch  angle  groups.  When  simulating  a  cold  beam,  this  number  should  be  one. 
Otherwise,  a  higher  number  is  needed.  Again,  the  number  needed  depends  on  many  factors  and  the 
best  way  to  judge  is  to  test  the  sensitivity  of  the  results  to  varying  the  number.  The  variable  ‘fntO’ 
controls  the  number  of  injected  particles  when  the  spectrum  has  multiple  frequencies.  If  there  are 
‘nw’  frequency  components  in  the  spectrum,  the  number  of  injected  particles  is  roughly 
‘na0*fnt0*nw*np0’.  For  a  single  frequency,  the  number  is  simply  ‘na0*np0’.  This  peculiar  way  of 
specifying  the  number  of  particles  is  designed  to  automatically  adjust  the  number  of  particles  to 
maintain  high  accuracy  as  the  parameters  of  the  RF  spectrum  are  changed.  The  setting  fntO  =1.5 
should  minimize  the  amount  of  aliasing  in  calculating  the  source  term  for  each  component  of  the 
spectrum.  However,  this  number  may  have  to  be  increased  if  the  spectrum  is  particularly  asymmetric 
about  ‘freqO’.  The  variable  ‘npO’  controls  the  basic  number  of  entrance  phases.  Its  size  depends  on 
the  number  of  space  charge  harmonics  that  are  believed  to  be  important.  At  a  minimum,  ‘npO’ 
should  be  on  the  order  of  20  (it’s  better  to  use  prime  numbers).  Also,  approximately  10  entrance 
phases  should  be  alloted  to  each  space  charge  harmonic.  The  number  of  space  charge  harmonics  is 
given  by  the  variable  ‘nsc’.  If  this  is  set  to  zero,  the  space  charge  field  is  not  included. 

3.2.4  Namelist  -  struct 

This  namelist  gives  some  general  parameters  of  the  structure.  Array  ‘zr’  is  an  ordered  set  of  axial 
positions  where  the  parameters  of  the  structure  are  to  be  specified.  The  first  element  must  be  zero  and 
the  last  must  be  greater  than  or  equal  to  the  length  of  the  structure,  which  is  given  by  variable  ‘zint’. 
Various  axially  dependent  parameters  will  be  defined  later  on  this  set  of  points.  Typically,  it  is 
assumed  that  quantities  vary  piecewise  linearly  between  points.  For  example,  ‘dvbeam’  is  a  vector  of 
voltages  defined  for  the  points  ‘zr’.  (Again,  the  first  element  should  be  zero.)  The  values  of 
‘dvbeam’  and  ‘zr’  determine  the  piecewise  linear  accelerating  potential  throughout  the  interaction 
region.  This  potential  can  act  to  accelerate  particles  to  improve  efficiency  similar  to  the  effect  of 
tapering  the  structure  period.  To  deactivate  this  potential,  set  all  elements  of  ‘dvbeam’  to  zero. 

The  number  of  entries  for  the  parameters  ‘zr’  and  ‘dvbeam’  should  not  exceed  by  more  than  one 
the  value  of  the  integer  parameter  ‘nzpts’  (this  parameter  is  declared  in  each  of  the  code’s 
subroutines).  That  is,  arrays  ‘zr’  and  ‘dvbeam’  are  declared  to  have  ‘nzpts+l’  elements.  The  same 
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will  be  true  for  a  number  of  variables  in  namelists  ‘helix_par’  and  ‘dis_dat’.  Applications  that  need 
many  axial  points  to  specify  the  profiles  of  the  structure  parameters  require  going  into  the  code  and 
changing  parameter  ‘nzpts’  in  each  place  that  it  appears.  When  fewer  than  ‘nzpts+l’  points  are 
required  to  specify  the  profiles,  it  is  usually  safe  to  input  fewer  than  that  many  values  for  the  relevant 
arrays.  This  is  illustrated  in  Section  4. 

Abrupt  changes  in  parameters  can  be  modeled  by  assigning  the  same  axial  position  to  two 
successive  elements  of  ‘zr’,  and  then  assigning  two  different  values  of  the  relevant  parameter  to  the 
appropriate  array.  The  code  will  place  the  values  at  two  adjacent  points  on  the  axial  grid.  It  is  not 
recommended  that  you  do  this  for  the  accelerating  potential  ‘dvbeam’  as  that  would  describe  an 
infinite  accelerating  field. 

The  logical  variable  ‘use_sh’  controls  the  way  the  dispersive  properties  of  the  structure  are  input 
to  the  code.  If  use_sh  =  .true.,  a  sheath  helix  model  is  assumed  and  the  parameters  for  the  sheath 
helix  are  given  in  namelist  ‘helix_par’.  If  use_sh  =  .false.,  the  impedance,  phase  velocity,  and 
equivalent  helix  radius  (required  to  calculate  the  space  charge  field)  are  read  from  namelist  ‘dis_dat’ 
or  read  in  from  files  ‘my .bp’  and  ‘my.zc’.  This  is  discussed  in  more  detail  in  the  example  described 
in  Section  4.8. 

3.2.5  Namelist  -  helix _par 

This  namelist  gives  arrays  of  the  parameters  of  the  sheath  helix;  ‘RH’  is  the  helix  radius  in  cm, 
‘H_LMDA’  is  the  helix  period  in  cm,  ‘RW’  is  the  wall  radius  in  cm,  ‘EPS_R’  describes  the  effective 
dielectric  constant  of  the  helix  supports,  and  ‘RV’  is  the  effective  vane  radius  in  cm.  If  there  are  no 
vanes,  ‘RV’  is  set  equal  to  wall  radius  ‘RW’.  The  various  values  are  defined  for  the  set  of  points  ‘zr’ 
given  in  namelist  ‘struct’.  The  number  of  entries  for  each  parameter  should  not  exceed  by  more  than 
one  the  value  of  the  integer  parameter  ‘nzpts’.  That  is,  the  code  expects  to  read  in  no  more  than 
‘nzpts+r  values  for  each  of  the  following  arrays:  ‘RH’,  ‘H_LMDA’,  ’RW’,  ‘EPS_R’,  and  ‘RV’. 
The  logical  ‘use_cor’  controls  the  form  of  the  space  charge  term  in  the  equation  of  motion.  If 
usejcor  -  .false.,  the  space  charge  term  is  calculated  with  the  assumption  that  the  helix  is  a  perfectly 
conducting  tube;  otherwise,  the  correction  to  the  space  charge  term  accounting  for  fields  outside  the 
helix,  as  calculated  in  Section  2.3,  is  included. 

3.2.6  Namelist  -  loss _par 

This  namelist  allows  one  to  specify  the  axial  dependence  of  the  attenuation  in  the  structure.  Here, 
‘za’  is  an  array  analogous  to  ‘zr’,  and  ‘att’  is  an  array  of  values  of  attenuation.  The  profile  of 
attenuation  is  taken  to  be  piecewise  linear  between  the  points  ‘za’.  The  maximum  allowed  length  of 
these  arrays  is  controlled  by  the  parameter  ‘nzapts’  defined  in  subroutines  ‘main’  and  ‘initparam’. 
The  values  of  attenuation  may  be  entered  either  in  dB/cm  or  in  dBA,  =  k^n  dB/(2;r)  (that  is,  dB  per 
axial  wavelength).  To  enter  the  data  in  dB/cm,  set  the  logical  ‘dB_cm’  to  true.  The  variable  ‘f_scale’ 
allows  for  flexibility  in  specifying  the  frequency  dependence  of  the  attenuation.  With  ‘f_scale’  =  0., 
the  frequency  dependence  of  the  attenuation  is  as  follows.  When  entering  in  dB/cm,  all  signal 
frequencies  will  have  the  same  rate  of  attenuation.  If  the  logical  ‘dB_cm’  is  set  to  false,  then  each 
signal  frequency  will  have  the  same  rate  of  attenuation  as  measured  in  dB/A,.  However,  since  different 
signal  frequencies  have  different  axial  wavelengths,  the  resulting  rates  of  attenuation  in  dB/cm  for 
each  signal  frequency  will  be  different.  With  ‘f_scale’  not  equal  to  0,  an  additional  frequency 
dependence  is  imposed  on  the  attenuation.  Specifically,  the  attenuation  for  each  frequency  /  is 
multiplied  by  a  factor  {f/freqOf-^’'^^ . 
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3.2.7  Namelist  -  ps _plot 

This  namelist  controls  some  of  the  output  from  CHRISTINE,  in  particular  information  about  the 
coordinates  of  particles  at  various  times  and  points  in  the  the  interaction  region.  Two  types  of  data 
files  for  making  plots  are  generated:  snapshots  of  the  positions  and  velocities  of  particles  at  given 
times,  and  phase  space  information,  namely  the  velocities  and  phases  of  particles  at  specified  axial 
positions.  These  files  are  not  written  if  the  program  is  in  the  parameter  scan  mode  {scan _p  =  .true.) 
The  variable  ‘nt_plot’  specifies  the  number  of  snapshot  files  to  be  written.  (If  nt _plot  =  0.,  then  no 
snapshot  files  are  written.)  As  the  fields  are  assumed  to  be  time  periodic,  snapshots  of  the  particle 
positions  and  velocities  are  time  periodic  as  well.  Thus,  it  does  not  make  sense  to  specify  an  absolute 
value  of  time  at  which  to  make  a  snap  shot.  Instead,  one  specifies  the  number  of  snapshots,  and  these 
are  taken  at  times  distributed  uniformly  over  the  period  of  the  fields.  If  only  one  frequency  is 
present,  ‘freqO’,  the  period  of  the  fields  is  the  reciprocal  of  this  frequency.  Otherwise,  the  period  is 
the  reciprocal  of  the  frequency  spacing  ‘delfreq’.  The  files  are  written  out  and  given  the  names 
‘bvsz_XXpYY’ ,  where  YY  is  the  number  of  snapshots  (‘nt_plot’)  and  XX  is  an  integer  between  1  and 
YY  indicating  the  time  interval.  The  first  of  these  plots  corresponds  to  the  reference  time  relative  to 
which  the  phases  of  the  input  signals  are  defined.  Equation  (1)  presents  the  convention  for  defining 
the  phases  of  the  signal  field. 

The  array  of  variables  ‘ps_z’  gives  the  axial  positions  at  which  phase  space  data  are  written.  If  a 
single  frequency  is  present,  the  phase  space  consists  of  the  particle  velocity  and  the  phase  y/Q.  If 
multiple  frequencies  are  present,  the  phase  space  is  the  particle  velocity  and  the  phase  Ay/.  The  phase 
space  files  are  given  the  names  ‘PS_ZZ’,  where  ZZ  is  the  index  corresponding  to  the  element  of 
array  ‘ps_z’. 

3.2.8  Namelist  -  dis_dat 

This  namelist  allows  one  to  enter  directly  the  relevant  parameters  for  simulating  a  structure  that  is 
not  well  described  by  a  sheath  helix,  or  one  for  which  measured  data  are  available.  The  variables 
‘npts’,  ‘f_in’,  ‘bp_in’,  and  ‘z_in’  determine  the  dispersive  properties  of  the  structure.  The  latter 
three  quantities  are  arrays  with  ‘ntpts’  elements.  Array  ‘f_in’  is  a  set  of  frequencies  that  spans  the 
range  ‘freq_min’  to  ‘freq_max’.  Array  ‘bp_in’  is  the  set  of  phase  velocities,  normalized  to  the 
speed  of  light,  corresponding  to  the  frequencies  in  ‘f_in’.  Array  ‘z_in’  is  the  corresponding  set  of 
coupling  impedances.  CHRISTINE  uses  these  arrays  to  interpolate  the  appropriate  values  of  phase 
velocity  and  impedance  for  the  set  of  frequencies  in  the  simulation. 

If  the  structure  is  tapered,  namelist  ‘dis_dat’  provides  a  quick  and  dirty  way  to  enter  the  tapering 
profile.  This  is  through  arrays  ‘fract_bp’  and  ‘fract_Z’.  These  arrays  allow  you  to  modify  the  data 
entered  in  arrays  ‘bp_in’,  and  ‘z_in’.  For  example,  ‘fract_bp’  is  an  array  of  factors  corresponding 
to  the  set  of  axial  points  defined  in  ‘zr’.  The  axially  dependent  phase  velocities  are  then  defined  by 
multiplying  the  value  of  ‘bp_in’  corresponding  to  the  frequency  of  interest  by  the  value  of 
‘fract_bp’  corresponding  to  the  axial  position  of  interest.  Similarly,  the  axially  dependent  impedance 
is  defined  as  the  product  of  ‘z_in’  and  ‘fract_Z’.  This  procedure  is  admittedly  crude  because 
varying  a  parameter  of  the  structure  does  not  result  in  the  phase  velocity  or  impedance  of  all 
frequencies  changing  by  the  same  factor.  However,  it  does  allow  for  a  treatment  of  tapered  profiles 
when  all  the  relevant  information  is  not  available. 

A  more  satisfactory  way  to  input  tapered  structure  data  is  to  set  the  logical  ‘use_data’  to  true.  The 
code  then  looks  for  files  named  ‘my .bp’  and  ‘my.zc’  containing  columns  of  numbers  describing  the 
phase  velocities  and  impedances  of  the  different  sections  of  the  structure.  In  this  case  the  frequency. 
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phase  velocity,  and  impedance  data  in  the  name  list  are  ignored.  The  format  of  the  files  ‘my.bp’  and 
‘my.zc’  should  be  as  follows.  There  should  be  ‘npts’  rows  of  data.  The  first  column  of  each  file 
should  contain  the  frequencies  replacing  ‘f_in’.  This  column  should  be  the  same  in  both  files.  The 
next  columns  in  each  file  contain  the  phase  velocities  or  impedances  for  each  of  the  points  defined 
by  ‘zr’  in  the  structure.  The  code  first  examines  ‘zr’  to  determine  how  many  columns  to  read.  It 
then  reads  in  the  data  via  an  unformatted  read  statement.  Thus,  tabs  or  spaces  may  be  used  to 
delineate  columns.  An  example  is  given  in  Section  4.8. 

3.2.9  Namelist  -  scan 

The  elements  of  this  namelist  control  the  parameter  scanning  feature  of  CHRISTINE.  To  activate 
the  scanning  feature,  set  ‘scan_p’  to  .true.  The  variables  that  can  be  scanned  are  listed  in  the 
comment  field  of  the  input  file.  CHRISTINE  will  then  perform  a  matrix  of  runs  corresponding  to  a 
scan  over  the  variables  ‘varX’  and  ‘varY’.  The  values  assigned  to  scanned  variables  then  override 
those  specified  in  the  other  namelists.  Additionally,  other  features  of  the  code  are  overridden  when  it 
is  in  the  scanning  mode.  First,  one  cannot  do  general  multifrequency  simulations  in  the  scanning 
mode.  Instead,  one  is  restricted  to  single  frequency,  or  multifrequency  with  harmonics-only 
simulations.  Specifically,  in  scanning  mode,  variables  ‘delfreq’  and  ‘freq_min’  are  set  equal  to  the 
relevant  value  of  ‘freqO’.  Variable  ‘freq_max’  is  set  equal  to  the  product  of  ‘freqO’  and 
‘n_harm’+l,  where  ‘n_harm’  is  the  number  of  harmonics  to  be  retained.  Thus,  to  simulate  a  single 
frequency,  set  ‘n_harm’  to  0.  The  values  that  the  variable  ‘varX’  will  take  on  are  controlled  by  the 
input  parameters  ‘nx’,  ‘varX_min’,  varX_max’,  and  ‘log_X’.  Basically,  ‘varX’  will  assume  ‘nx’ 
values  between  ‘varX_min’  and  ‘varX_max’,  inclusive.  If  ‘nx’  =  1,  only  the  minimum  value  is 
simulated.  If  ‘log_X’  is  .false.,  the  values  of  ‘varX’  will  be  incremented  linearly  between  the  two 
limits.  With  ‘log_X’  set  to  .true.,  the  logarithm  of  ‘varX’  is  incremented  linearly.  Similar  rules  apply 
to  ‘varY’. 

In  the  scanning  mode,  CHRISTINE  produces  the  following  output  files:  ‘gain’,  ‘efficiency’,  and 
‘P_out.N’.  (Files  such  as  ‘axial’  and  those  containing  particle  trajectory  data  are  not  printed  out 
when  the  code  is  in  the  parameter  scanning  mode.)  Each  of  the  output  files  contains  a  table  where  the 
rows  correspond  to  the  different  values  of  ‘VarX’  and  the  columns  correspond  to  the  different  values 
of  ‘varY’.  The  file  ‘gain’  contains  the  magnitude  and  phase  of  the  gain  of  the  fundamental 
frequency  signal  for  each  of  the  input  parameters.  The  files  ‘P_out.N’  contain  the  output  power  in 
the  ‘N  th’  signal  frequency,  where  N  =  1  is  the  fundamental,  N  =  2  the  harmonic,  and  so  on.  The  file 
‘efficiency’  contains  the  efficiency  of  beam-to-output  power  conversion  based  on  the  output  power 
in  the  fundamental  frequency. 

4.  SAMPLE  RUNS 

This  section  presents  several  illustrations  of  the  code’s  features  by  providing  sample  namelist  files. 
The  physical  parameters  used  are  based  on  a  specific  research  tube.  The  following  examples  are 
presented: 

1.  Sheath  helix  modeling 

2.  Phase  space  plots 

3.  Scanning  parameters 

4.  Space  charge 

5.  Harmonic  generation 

6.  Multifrequency  simulations 

7.  Tapering  and  attenuation 

8.  Importing  dispersion  and  coupling  impedance  data. 
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The  first  several  examples  use  the  input  file  shown  in  Table  1 . 

4.1  Sheath  Helix  Modeling 

The  sheath  helix  model  is  activated  by  setting  the  logical  variable  ‘use_sh’  to  true  in  namelist 
‘struct’.  The  parameters  of  the  helix  are  then  entered  in  namelist  ‘helix_par’.  (If  ‘use_sh’  is  set  to 
false,  then  none  of  the  information  in  ‘helix_par’  is  used.  Instead,  information  in  ‘dis_dat’  is  used. 
This  is  discussed  in  Section  4.8.)  For  this  example,  an  untapered  helix  with  the  following  dimensions 
is  specified:  helix  radius  =  0.12446  cm,  helix  period  =  0.080137  cm,  wall  radius  =  0.2794  cm,  vane 
radius  =  0.2794  cm,  and  effective  dielectric  constant  of  helix  supports  =  1.75.  Note  that  in  this  case, 
the  wall  radius  and  vane  radius  are  the  same.  This  corresponds  to  the  case  with  no  vanes.  Note  also 
that  two  entries  appear  for  each  parameter  (e.g.,  RH  =  2*.12446).  This  is  because  the  axial  profile  of 
the  corresponding  variable  is  taken  to  be  a  piecewise  linear  function  between  successive  values  of  the 
elements  of  array  ‘zr’  appearing  in  namelist  ‘struct’.  In  the  current  example,  the  first  two  elements  of 
‘zr’  are  0.0  and  10.0.  As  this  second  number  is  the  first  to  exceed  the  length  of  the  interaction  region 
{zint  =  9.5758),  all  subsequent  elements  of  ‘zr’  are  unimportant  and  the  axial  profile  of  various 
functions  is  taken  to  be  piecewise  linear  for  z  between  0.0  cm  and  10.0  cm.  In  order  to  represent  a 
constant,  one  must  specify  the  two  values  of  the  parameters  corresponding  to  the  points  zr  =  0  and 
zr  =  10  to  be  the  same. 

The  characteristics  of  the  RF  spectrum  are  controlled  by  namelist  ‘rad’.  The  present  case 
prescribes  a  5  GHz,  single  frequency  run.  The  minimum  and  maximum  frequencies  (‘freq_min’  and 
‘freq_max’)  are  the  same  as  the  central  frequency  (‘freqO’)  and  the  frequency  separation 
(‘delfreq’).  Examples  of  multifrequency  runs  are  given  in  Section  4.6.  The  input  power  of  the  signal 
(‘pow_in’)  is  2.4  x  lO’^W,  and  the  phase  (phase_in)  is  zero  radians.  As  discussed  in  Section  2.2,  the 
phase  is  defined  with  respect  to  a  reference  time  to  which  the  particle  diagnostics  discussed  in  Section 

4.2  are  keyed.  The  ballistic  prebunching  amplitude  and  phase  are  both  zero,  and  the  current  is 
ungated  (i.e.,  gated_l  =  .false.). 

Numerical  parameters  are  controlled  by  namelist  ‘num’  (Section  3.2.3).  Sheath  helix  modeling 
offers  a  total  of  137  entrance  times,  which  is  probably  many  more  than  is  necessary.  Only  one  space 
charge  harmonic  is  included  in  the  simulation.  The  effect  of  varying  the  number  of  space  charge 
harmonics  is  discussed  in  Section  4.4. 

Running  the  code  produces  the  following  output  files:  ‘dispersion’,  ‘impedance’,  ‘axial’, 
‘P_out’,  ‘bvsz_01p03’,  ‘bvsz_02p03’,  ‘bvsz_03p03’,  ‘PS.Ol’,  ‘PS.02’,  and  ‘PS.03’.  Files 
‘dispersion’  and  ‘impedance’  contain  tables  of  the  normalized  phase  velocity  and  coupling 
impedance,  respectively,  calculated  in  the  sheath  helix  approximation  for  the  frequency  spectrum 
prescribed  in  namelist  ‘rad’.  This  example  offers  only  one  frequency  (5  GHz).  CHRISTINE 
tabulates  two  phase  velocities  and  coupling  impedances  for  each  frequency,  corresponding  to  the 
values  of  the  relevant  quantity  at  the  set  of  axial  points  contained  in  ‘zr’.  In  this  example,  therefore, 
two  identical  values  are  tabulated  for  each  quantity.  The  example  described  in  Section  4.6  generates  a 
more  interesting  table. 

File  ‘axial’  contains  a  table  of  the  axial  profiles  of  the  average  beam  energy  and  the  power  in  each 
frequency  component  of  the  RF  signal.  Figure  2  plots  the  data  for  the  present  run.  File  ‘P_out’ 
contains  a  table  of  output  power  vs  frequency.  This  information  is  already  present  in  the  last  line  of 
the  file  ‘axial’,  and  thus,  is  redundant.  However,  when  doing  multifrequency  simulations,  use  of  this 
table  provides  a  convenient  way  of  displaying  the  output  spectrum. 
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Fig.  2  —  File  ‘axial’:  power  and  average  beam  energy  vs  z 


Files  ‘bvsz_01p03’.  ‘bvsz_02p03’,  ‘bvsz_03p03’,  ‘PS.Ol’,  ‘PS.02’.  and  ‘PS.03’  display 
information  about  the  trajectories  of  individual  particles  and  are  discussed  in  Section  4.2. 

4.2  Phase  Space  Plots 

Namelist  ‘ps_plot’  controls  the  output  of  the  files  containing  data  on  the  particle  trajectories.  The 
parameter  ‘nt_plot’  gives  the  number  of  files  containing  the  axial  position  and  normalized  axial 
velocity  coordinates  of  particles  at  a  fixed  time.  This  example  produces  three  such  files: 
‘bvsz_01p03’,  ‘bvsz_02p03’,  and  ‘bvsz_03p03’.  The  first  of  these  gives  the  coordinates  at  a  time 
corresponding  to  the  reference  time  (keyed  to  the  phase  of  the  input  signal)  and  the  second  and  third 
give  the  coordinates  of  the  particles  one  third  and  two  thirds  of  a  period  of  the  signal  later, 
respectively.  Figure  3  is  a  scatter  plot  of  the  data  in  ‘bvsz_01p03’. 

Files  ‘PS.Or,  ‘PS.02’,  and  ‘PS.03’  contain  phase  space  data  (i.e.,  tables  of  phase  and  normalized 
axial  velocity  of  particles  at  specific  values  of  axial  position).  The  axial  positions  are  controlled  by 
the  elements  of  array  ‘ps_z’.  The  files  are  numbered  consecutively  corresponding  to  these  elements. 
Figure  4,  a  plot  of  the  phase  space  coordinates  contained  in  ‘PS.03’  corresponding  to  z  =  9  cm, 
clearly  shows  the  saturation  of  the  amplifier  caused  by  phase  trapping. 

Particle  trajectory  data  are  not  printed  out  when  the  code  is  in  the  parameter  scanning  mode,  which 
is  illustrated  in  Section  4.3. 
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4.3  Scanning  Parameters 

The  preceding  example  produced  data  for  a  single  set  of  input  parameters.  Often  it  is  desired  to 
vary  parameters  and  plot  the  results  vs  the  parameter  that  is  varied.  This  feature  is  controlled  in 
CHRISTINE  by  namelist  ‘scan’.  To  activate  the  parameter  scanning  feature,  the  logical  parameter 
‘scan_p’  is  set  to  true.  A  sample  scan  namelist  is  shown  below. 

$scan 

scan_p  =  .tme. 
n_harm  =  0, 
varX  =  ‘vbeam’, 
nx  =  50, 

varX_min  =  2.0e3, 
varX_max  =  4.0e3, 
log_X  =  .false. 
varY  =  ‘freqO’, 
ny  =  2, 

varY_min  =  5.e9, 
varY_max  =  6.e9, 
log_Y  =  .false. 

Send 


This  namelist  prescribes  that  two  variables,  ‘vbeam’  and  ‘freqO’,  are  to  be  varied.  There  will  be  50 
different  values  of  ‘vbeam’  between  2.0  kV  and  3.0  kV,  and  two  values  of  ‘freqO’  between,  in  this 
example,  5  GHz  and  6  GHz.  The  code  will  thus  simulate  the  operation  of  the  device  for  100  different 
sets  of  parameters. 

Running  the  code  with  the  sample  namelists  produces  the  following  output  files:  ‘gain’, 
‘efficiency’,  and  ‘P_out.r.  (Files  such  as  ‘axial’  and  those  containing  particle  trajectory  data  are 
not  printed  out  when  the  code  is  in  the  parameter  scanning  mode.)  Each  of  these  contains  a  table 
where  the  rows  correspond  to  the  different  values  of  ‘VarX’  (in  this  case  ‘vbeam’)  and  the  columns 
correspond  to  the  different  values  of  ‘VarY’  (in  this  case  ‘freqO’).  File  ‘gain’  contains  the 
magnitude  and  phase  of  the  gain  of  the  fundamental  frequency  signal  for  each  of  the  input 
parameters.  File  ‘P_out.r  contains  the  output  power  in  the  fundamental  frequency,  and  ‘efficiency’ 
contains  the  efficiency  of  the  beam-to-energy  power  conversion  based  on  the  output  power  in  the 
fundamental  frequency.  If  harmonics  are  included  by  setting  ‘n_harm’  to  an  integer  different  from 
zero,  then  the  output  power  in  each  harmonic  is  printed  in  file  ‘P_out.N’,  where  N  is  an  integer 
labeling  the  frequency  of  the  corresponding  signal,  1  implies  the  fundamental,  2  the  second 
harmonic,  and  so  on.  Figure  5  plots  file  ‘P_out.r. 
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Fig.  5  —  File  ‘P_out.r:  output  power  vs  beam  voltage  for  two  values  of  the  input  frequency 


4.4  Space  Charge 

The  number  of  space  charge  harmonics  is  controlled  by  parameter  ‘nsc’  in  namelist  ‘num’ .  As 
discussed  in  Section  4.1,  the  space  charge  field  can  have  a  strong  temporal  harmonic  content  even  if 
the  RF  signal  field  is  relatively  pure.  Figure  6  illustrates  this  by  displaying  the  results  of  a  number  of 
parameter  scan  runs  with  differing  numbers  of  space  charge  harmonics.  (Parameter  ‘nsc’  cannot  be 
scanned  automatically,  so  this  plot  is  made  by  collecting  data  from  a  number  of  runs.)  Generally 
speaking,  as  the  number  of  space  charge  harmonics  increases,  so  does  the  number  of  particles  ‘npO’. 
In  this  example,  all  runs  were  conducted  with  ‘npO’  =  137.  Additionally,  the  sheath  helix  correction 
to  the  space  charge  field  was  not  activated  (i.e.,  use_cor  =  .false.). 


Fig.  6  —  Output  power  vs  beam  voltage  for  ‘freqO’  =  5.0  GHz  and  several  values  of  ‘nsc’ 
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The  plot  shows  that  between  five  and  ten  space  charge  harmonics  are  required  to  obtain  a  converged 
result.  This  is  because  the  output  has  saturated  and  the  bunched  beam  current  has  a  high  harmonic 
content. 

4.5  Harmonic  Generation 

This  example  illustrates  a  multifrequency  simulation  for  a  single  set  of  parameters  (i.e.,  scan^  = 
.false.).  The  spectrum  is  specified  by  setting  the  value  of  ‘freq_max’  to  20.e09  (i.e.,  20  GHz).  With 
‘freqO’,  ‘delfreq’,  and  ‘freq_min’  all  equal  to  5.e09,  the  spectrum  will  consist  of  the  following 
frequencies;  5  GHz,  10  GHz,  15  GHz,  and  20  GHz.  Namelist  Tad’  appears  below. 

$rad 

freqO  =  5.e9, 
delfreq  =  5.e9, 
freq_max  =  20.e9, 
freq_min  =  5.e9, 
pow_in  =  2.4e-03,3*0. 
phase_in  =  4*0., 

QB  =  4*.0, 
phase_B  =  4*0., 
gated_I  =  .false., 
sin_2  =  .true., 

Avel_pkl  =  .60, 

Send 

Note  that  we  have  increased  the  number  of  entries  for  arrays  ‘powjn’,  ‘phase_in’,  ‘QB’,  and 
‘Phase_B’  so  that  it  corresponds  to  the  number  of  signals  in  the  spectrum.  Had  we  desired  to  inject  a 
signal  at  the  second  harmonic  (10  GHz),  we  would  have  entered  the  corresponding  values  of  input 
power  and  phase  as  the  second  elements  of  arrays  ‘pow_in  and  phase_in  ,  respectively. 
Additionally,  for  this  run  we  included  five  space  charge  harmonics.  Figure  7  plots  ‘axial’. 
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Fig.  7  —  File  ‘axial’:  power  in  various  frequencies  and  average  beam  energy  vs  z 


4.6  Multifrequency  Simulation 

This  example  shows  the  results  of  a  large  scale  multifrequency  simulation;  CHRISTINE  was 
designed  for  just  this  type  of  simulation.  Namelist  ‘rad’  is  modified  as  shown  below. 

$rad 

freqO  =  5.e9, 
delfreq  =  l.e9, 
freq_max  =  15.e9, 
freq_min  =  1  .e9, 
powjn  =  4*0.,2*2.4e-03,9*0. 
phase_in  =  15*0., 

QB  =15*.0, 
phase_B  =  15*0., 
gated_I  =  .false., 
sin_2  =  .true., 

Avel_pkl  =  .60, 

Send 


The  spectrum  consists  of  15  frequencies  spaced  1  GHz  apart  starting  from  1  GHz  and  extending 
to  15  GHz.  Thus,  15  entries  are  required  for  each  array  connected  with  the  input  signal.  In  this  case, 
the  fifth  and  sixth  elements  of  the  input  power  corresponding  to  5  GHz  and  6  GHz  are  driven  at  2.4  x 
10'^  W.  Not  shown  is  namelist  ‘num’,  where  parameter  ‘npO’  has  been  reduced  to  71.  With  the  15 
frequency  components,  a  total  of  1562  entrance  times  are  available.  If  more  entrance  times  are 
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required,  the  value  of  parameter  ‘nnp’  that  appears  in  the  parameter  statement  of  each  subroutine 
must  be  increased. 

Figure  8  plots  the  contents  of  file  ‘P_out’.  The  intermodulation  between  the  5  GHz  and  6  GHz 
signal  is  apparent  in  that  signals  are  present  at  every  integer  frequency. 


frequency  [GHz] 

Fig.  8  —  Contents  of  output  spectrum  ‘P_out’ 


Figure  9  is  an  example  of  a  phase  space  plot  (‘PS.03’)  for  this  multifrequency  simulation.  Here, 
variable  6,  which  spans  In,  represents  the  arrival  time  of  a  particle  at  z  =  9.0  cm  (ps_z(3)  =  9.0)  times 
the  angular  frequency  separation  Am.  The  dominance  of  the  5  GHz  signal  in  the  output  (Fig.  8)  is 
evidenced  by  the  apparent  m  =  5  periodicity  in  the  phase  space. 


-4  -3-2-1  0  1  2  3  4 


theta 

Fig.  9  —  Contents  of  ‘PS.03’  showing  phase  space  at  z  =  9.0  cm 
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Finally,  Fig.  10  displays  on  a  single  set  of  axes  the  content  of  files  ‘dispersion’  and  ‘impedance’. 
These  show  the  characteristics  of  the  structure  as  calculated  in  the  sheath  helix  approximation  (see 
Section  4.1). 


Fig.  10 —  Phase  velocity  and  impedance  vs  frequency 


4.7  Tapering 

An  example  of  a  tapered  structure  with  attenuation  is  illustrated  by  making  the  following  changes 
to  namelists  ‘struct’,  ‘helix_par’,  and  ‘loss_par’. 

Sstruct 

zint  =  9.5758, 
zr  =  0.  ,7.,  8.,  10., 
dvbeam  =  4*0. 
use_sh  =  .true., 

Send 

$helix_par 

RH  =  4*.12446, 

H_LMDA  =  2*.080137,.07212,.06491 
RW  =  4*.2794, 

EPS_R  =  4*1.75, 

RV  =  4*.2794, 
use_cor  =  .false. 

Send 
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$loss_par 

za  =  0.,2.,4.,6.,8.,10., 
att  =2*1.,5.,3*1. 
dB_cm  =  .true. 
f_scale  =  0. 

Send 

The  parameters  of  the  structure  are  defined  at  four  different  points:  z  =  0,  7,  8,  and  10  cm.  The 
parameter  that  varies  is  the  helix  period  ‘H_LMDA’.  Between  z  =  0  and  z  =  7  cm,  it  is  equal  to 
0.080137  cm.  It  then  varies  linearly  to  a  value  of  0.07212  cm  at  z  =  8  cm  and  then  to  0.06491  at  z  = 
10  cm.  Attenuation  is  similarly  defined.  The  attenuation  is  1  dB/cm  for  0  <  z  <  2  cm.  It  then  makes  a 
linear  tansition  to  a  value  of  5  db/cm  at  z  =  4  cm.  It  then  transitions  back  to  1  dB/cm  by  z  =  6  cm  and 
remains  at  that  value  throughout  the  remainder  of  the  interaction  region. 

Files  ‘dispersion’  and  ‘impedance’  now  contain  a  number  of  columns  giving  the  frequency 
dependence  of  the  relevant  quantities  at  each  of  the  four  points  (z  =  0,  7,  8,  and  10  cm).  Table  2 
shows  the  output  of  file  ‘dispersion^ 

Due  to  the  attenuation,  the  device  output  power  is  somewhat  lower  than  in  the  previous  examples. 
Figure  1 1  plots  the  contents  of  file  ‘P_out’ . 


Table  2 —  Contents  of  File  ‘Dispersion’:  Phase  Velocity  vs  Frequency  at  Four  Points 

Along  the  Structure 


freq  [GHz] 

0 

1 

2 

3 

.lOOOE+01 

.1079E+00 

.1079E+00 

.9714E-01 

.8735E-01 

.2000E+01 

.1053E+00 

.1053E+00 

.9428E-01 

.8425E-01 

.3000E+01 

.1016E+00 

.1016E+00 

.9045E-01 

.8034E-01 

.4000E+01 

.9779E-01 

.9779E-01 

.8672E-01 

.7687E-01 

.5000E+01 

.9448E-01 

.9448E-01 

.8381E-01 

.7441E-01 

.6000E+01 

.9200E-01 

.9200E-01 

.8179E-01 

.7284E-01 

.7000E+01 

.9027E-01 

.9027E-01 

.8047E-01 

.7187E-01 

.8000E+01 

.8909E-01 

.8909E-01 

.7962E-01 

.7128E-01 

.9000E+01 

.8830E-01 

.8830E-01 

.7907E-01 

.7090E-01 

.lOOOE+02 

.8777E-01 

.8777E-01 

.7871E-01 

.7067E-01 

.llOOE+02 

.8740E-01 

.8740E-01 

.7847E-01 

.7051E-01 

.1200E+02 

.8714E-01 

.8714E-01 

.7831E-01 

.7041E-01 

.1300E+02 

.8696E-01 

.8696E-01 

.7819E-01 

.7034E-01 

.1400E+02 

.8684E-01 

.8684E-01 

.7812E-01 

.7030E-01 

.1500E+02 

.8675E-01 

.8675E-01 

.7806E-01 

.7026E-01 

Fig.  11  —  File  ‘P_out’:  the  output  spectrum  for  the  case  with  tapering  and  attenuation 


4.8  Importing  Data 

Instead  of  using  the  sheath  helix  model  to  calculate  the  dispersive  properties  of  the  structure,  one 
may  wish  to  import  data  for  a  particular  structure.  This  is  done  by  setting  the  logical  variable 
‘use_sh’  in  ‘struct’  to  ‘false’.  Data  can  now  be  entered  through  namelist  ‘dis_dat’,  or  if  the 
dispersion  data  is  contained  in  files  (named  ‘my .bp’  and  ‘my.zc’),  these  may  be  read  in  as  well. 

We  first  discuss  the  case  in  which  the  relevant  data  is  entered  directly  into  namelist  ‘dis_dat’.  A 
sample  appears  below.  Data  is  entered  in  the  form  of  a  set  of  arrays  giving  frequencies,  phase 
velocities,  and  impedances,  each  with  ‘npts’  entries.  The  frequencies  need  not  correspond  one  to  one 
to  the  input  signals  specified  in  ‘rad’,  however  they  must  span  the  range  of  frequencies  ‘freq_min’ 
to  ‘freq_max’.  The  code  will  interpolate  phase  velocity  and  impedance  between  given  values  if 
necessary.  If  the  structure  is  tapered,  namelist  ‘dis_dat’  enables  the  user  to  enter  the  tapering  profile" 
through  the  use  of  arrays  ‘fract.bp’  and  ‘fract_Z’.  Below  is  an  example  of  a  namelist  using  the 
dispersion  and  impedance  data  displayed  in  Fig.  11  along  with  a  quick  and  dirty  (Q&D)  prescription 

of  tapering  designed  to  model  the  previous  example  that  illustrated  tapering'  using  the  sheath  helix 
model. 

$dis_dat 
npts  =  15, 
f_in  = 

I.e09,2.e09,3.e09,4.e09,5.e09,6.e09,7.e09,8.e09,9.e09,10.e09,ll.e09,12.e09,13.e09,14e09,15.e09, 

bp_in  = 

0.10790,0.10530,0.10160,0.097790,0.094480,0.092000,0.090270,0.089090,0.088300,0.087770,0.0 

87400,0.087140,0.086960,0.086840,0.086750 

z_in  = 
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236.40,212.30,172.00,123.20,78.920,47.030,27.090,15.430,8.8040,5.0490,2.9170,1.6970,0.99390,0 

.58570,0.34710 

fract_bp  =2*1.,.89,.79 
fract_Z  =  2*1.,.84,.68 
RH_eff  =  .12446, 
use_data  =  .false., 

$end 

The  values  of  array  ‘fract_bp’  are  chosen  so  that  the  profile  of  the  phase  velocity  at  5  GHz 
matches  that  displayed  in  the  table  in  the  previous  example.  Similarly,  the  elements  of  ‘fract_Z’  are 
picked  to  match  the  coupling  impedance  at  5  GHz  of  the  previous  example.  Figure  12  provides  a 
comparison  of  the  dispersion  and  coupling  impedance  at  the  end  of  the  structure  using  the  sheath 
helix  model  with  tapered  parameters  (the  previous  example)  and  the  imported  data  with  the  Q&D 
tapering  procedure. 
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Fig.  12 —  Comparison  of  structure  data  for  sheath  helix  model  and  Q&D  tapering 


As  can  be  seen,  the  Q&D  method  succeeds  as  designed  in  producing  the  correct  values  of 
impedance  and  phase  velocity  at  5  GHz.  However,  the  parameters  do  not  match  at  other  frequencies. 
The  effect  this  has  on  the  output  spectrum  is  displayed  in  Fig.  13. 

If  the  Q&D  method  does  not  meet  the  user’s  needs,  tapered  dispersion  and  impedance  data  may 
be  imported  directly.  To  do  this,  the  logical  ‘use_data’  must  be  set  to  true  and  the  code  must  be 
supplied  with  the  two  files  ‘my .bp’  and  ‘my.zc’  in  the  format  shown  in  Table  3.  As  an  example,  the 
table  can  be  used  to  import  the  ‘correct’  sheath  helix  data  into  the  code;  the  user,  of  course,  would 
supply  data  generated  by  some  other  means. 
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Fig.  13  —  Comparison  of  output  spectrum  for  two  methods  of  introducing  tapered  parameters 

Table  3  —  Sheath  Helix-generated  Data  in  the  Format  to  be  Read  as  ‘my.bp’ 

freq  [GHz]  0  1  2  3 

.lOOOE+10  .1079E+00  .1079E+00  .9714E-01  .8735E-01 

.2000E+10  .1053E+00  .1053E+00  .9428E-01  .8425E-01 

.3000E+10  .1016E+00  .1016E+00  .9045E-01  .8034E-01 

.4000E+10  .9779E-01  .9779E-01  .8672E-01  .7687E-01 

.5000E+10  .9448E-01  .9448E-01  .8381E-01  .7441E-01 

.6000E+10  .9200E-01  .9200E-01  .8179E-01  .7284E-01 

.7000E+10  .9027E-01  .9027E-01  .8047E-01  .7187E-01 

.8000E+10  .8909E-01  .8909E-01  .7962E-01  .7128E-01 

.9000E+10  .8830E-01  .8830E-01  .7907E-01  .7090E-01 

•  lOOOE+ll  .8777E-01  .8777E-01  .7871E-01  .7067Er01 

•  llOOE+ll  .8740E-01  .8740E-01  .7847E-01  .7051E-01 

.1200E+11  .8714E-01  .8714E-01  .7831E-01  .7041E-01 

.1300E+11  .8696E-01  .8696E-01  .7819E-01  .7034E-01 

.1400E+11  .8684E-01  .8684E-01  .7812E-01  .7030E-01 

.1500E+11  .8675E-01  .8675E-01  .7806E-01  .7026E-01 


The  elements  of  Table  3  (‘my.bp’)  are  separated  by  blank  spaces  and  each  line  is  terminated  with 
a  carriage  return.  The  code  then  produces  results  that  are  identical  to  those  of  the  tapering  example 
using  the  sheath  helix  option. 
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